This file and all other referenced in the code can be found at the repository: https://github.com/Andros-Spica/ERL_Angourakis-et-al-2020_Rproject
Using:
Daniel Wallach, David Makowski, James W. Jones, François Brun (2019), ‘Working with dynamic crop models: Methods, tools, and examples for agriculture and enviromnent’, Academic Press.
Model description in p. 24-28, R code example in p. 116-122. Another, possibly newer version of the R code can be found at: https://github.com/cran/ZeBook/blob/master/R/watbal.model.r
Original model from:
Woli P, Jones J W, Ingram K T and Fraisse C W (2012) ‘Agricultural Reference Index for Drought (ARID)’ Agron. J. 104-287. Online: https://www.agronomy.org/publications/aj/abstracts/104/2/287
Load functions from ‘Soil Water Balance model’ library file:
source("source/watbal.model.R")
Technical note: The ZeBook package, watbal.model.R has a bug in lines 21-25 and 58. The problem is that the code fails to extract values from the param object generated by watbal.define.param because it does not account for both dimensions of this object: (WHC, MUF, DC, z, CN) x (nominal, binf, bsup). I solve this problem by adding “typeOfParameterValue” as an argument in watbal.model that is passed also to watbal.update.
Load function to estimate reference evapotranspiration:
source("source/estimateETr.R")
Load CurlyBraces function for plotting:
source("source/CurlyBraces.R")
< begin parentesis >
To reproduce Figure 4.5 in Wallach et al. 2006, load weather daily data and select site and year :
weather <- watbal.weather(working.year = 2008, working.site = 1)
Or following code on the book:
weather <- ZeBook::weather_FranceWest
weather <- weather[(weather$idsite == 1 & weather$WEYR == 2008),]
And skip to the steps related to the Soil Water Balance model (bellow).
< end parentesis >
To generate Figure 3, we use the data downloaded at NASA´s POWER access viewer (power.larc.nasa.gov/data-access-viewer/) selecting the user community ‘Agroclimatology’ and pin pointing the five different locations between 01/01/1984 and 31/12/2018. The exact locations are:
We selected the ICASA Format’s parameters:
The following code reads all files inside the “Fig3_input” folder and creates a single data frame that contains the data of all five locations:
inputFiles <- paste0("Fig3_input/", list.files(path = "Fig3_input"))
weather <- data.frame()
for (i in 1:length(inputFiles))
{
tempData <- watbal.weather.file(inputFiles[i], year = 1984:2008)
# convert any missing data in all sky insolation (I) from -99.0 to Na
tempData$I[tempData$I == -99] <- NA
weather <- rbind(weather, tempData)
}
rm(tempData)
# use Latitude to distinguish sites and create a new variable using site names:
weather$Site <- rep(NA, nrow(weather))
for (i in 1:nrow(weather))
{
if (floor(weather$LAT[i]) == 27) { weather$Site[i] <- "Mohenjo-daro" }
if (floor(weather$LAT[i]) == 23) { weather$Site[i] <- "Dholavira" }
if (floor(weather$LAT[i]) == 28) { weather$Site[i] <- "Ganweriwala" }
if (floor(weather$LAT[i]) == 30) { weather$Site[i] <- "Harappa" }
if (floor(weather$LAT[i]) == 29) { weather$Site[i] <- "Rakhigarhi" }
}
weather$Site <- factor(weather$Site,
levels = c("Mohenjo-daro", "Ganweriwala", "Harappa",
"Dholavira", "Rakhigarhi"))
Sum up yearly and seasonal totals by splitting the year into summer (May-Oct or from day ) and winter (Nov-Apr):
precipitationSums <- list(Site = c(), year = c(), total = c(), summer = c(), winter = c())
for (site in levels(weather$Site))
{
for (year in levels(factor(weather$year)))
{
precipitationSums$Site <- c(precipitationSums$Site, site)
precipitationSums$year <- c(precipitationSums$year, year)
precipitationSums$total <- c(precipitationSums$total,
round(sum(weather$RAIN[weather$Site == site &
weather$year == year])))
precipitationSums$winter <- c(precipitationSums$winter,
round(sum(weather$RAIN[weather$Site == site &
weather$year == year &
(weather$day <= 121 | weather$day > 305)]
)))
precipitationSums$summer <- c(precipitationSums$summer,
round(sum(weather$RAIN[weather$Site == site &
weather$year == year &
(weather$day > 121 & weather$day <= 305)]
)))
}
}
precipitationSums <- data.frame(precipitationSums)
Calculate the average annual precipitation in every site to reference in the text:
for (site in levels(precipitationSums$Site))
{
cat(paste0(site, "\n"))
print(summary(precipitationSums$total[precipitationSums$Site == site], na.rm = TRUE))
cat("---------------------------------------\n")
}
## Dholavira
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 55.0 188.0 258.0 303.3 420.0 671.0
## ---------------------------------------
## Ganweriwala
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 30 85 107 117 149 219
## ---------------------------------------
## Harappa
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 123 164 232 240 299 459
## ---------------------------------------
## Mohenjo-daro
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 8.00 33.00 72.00 78.96 98.00 286.00
## ---------------------------------------
## Rakhigarhi
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 248 392 495 499 615 740
## ---------------------------------------
Estimate reference evapotranspiration using FAO recomendations (Penman-Monteith method):
weather$ETr <- estimateETr(weather$I,
weather$T2M,
weather$Tmax,
weather$Tmin,
weather$T2MDEW,
weather$WS2M)
Use Soil Water Balance model to calculate soil water content proportion or WATp (mm mm-1 day-1) and ARID coefficient.
The first step is to set all soil parameters, here assumed to be constants throughout the sites:
parameters <- watbal.define.param()
# This function creates the following configuration (from Wallach et al. 2006):
# WHC = 0.15
# MUF = 0.096
# DC = 0.55
# z = 40
# CN = 65
# input variables describing soil (from Wallach et al. 2006)
soil.WP = 0.06
soil.FC = soil.WP + parameters["nominal", "WHC"]
Run the model inputting parameters and dataset:
aWatbal <- cbind(Site = weather$Site,
watbal.model(parameters,
weather,
WP = soil.WP,
FC = soil.FC,
typeOfParValue = "nominal"))
watbal.modeliterates for every day in the dataset calculating soil water content, absolute and proportion (WAT, WATp), and the respective ARID coefficient.
Plot precipitation (RAIN), evapotranspiration (ETr), soil water content proportion (WATp), and the ARID coefficient for the five sites during the year 1990.
selectedYear = 1990
lastDayOfWinter = 121
lastDayOfSummer = 305
#--------
# To create a png file:
plotName = "Fig3.png"
fontScale = 1.2
png(plotName, width = 1000, height = 600)
#---------
# alternatively, to create eps file:
# plotName = "Fig3.eps"
#
# fontScale = 1
#
# extrafont::loadfonts(device = "postscript")
# grDevices::postscript(file = plotName,
# pointsize = 10,
# width = 10,
# height = 6,
# horizontal = FALSE,
# paper = "special",
# onefile = FALSE,
# family = "sans",
# colormodel = "cmyk")
#---------
layout(matrix(1:25, nrow = 5, ncol = 5),
heights = c(0.2, 1.1, 1, 1, 1.3),
widths = c(1.25, 1, 1, 1, 1))
yLabs <- c("Precipitation (mm)", "ETr (mm)", "WATp (%)", "ARID index")
leftMargin <- 4.1
# x coordinate in barplots has different scale, so must be adjusted with this offset
curlyBracesOffset = 1.2
for (site in levels(factor(weather$Site)))
{
if (site != levels(factor(weather$Site))[1])
{
yLabs <- rep("", 4)
leftMargin <- 2.1
}
tempData <- aWatbal[aWatbal$year == selectedYear & aWatbal$Site == site,]
par(mar = c(0, leftMargin, 0, 0), cex = fontScale)
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
text(x = 0.5, y = 0.7, font = 4, cex = 1.2,
labels = site)
text(x = 0.5, y = 0.1, cex = 0.8,
labels = paste0(precipitationSums$total[precipitationSums$Site == site &
precipitationSums$year == selectedYear],
" mm (annual total)")
)
par(mar = c(1.1, leftMargin, 2.1, 0.2))
barplot(tempData$RAIN,
ylab = yLabs[1])
CurlyBraces(x0 = (lastDayOfWinter + 1) * curlyBracesOffset,
x1 = lastDayOfSummer * curlyBracesOffset,
y0 = 1.1 * max(tempData$RAIN), y1 = 1.1 * max(tempData$RAIN),
position = 3,
bracesSize = 0.15 * max(tempData$RAIN),
label = paste0(
precipitationSums$summer[precipitationSums$Site == site &
precipitationSums$year == selectedYear],
" mm (summer)"),
labelSize = 1, labelDist = 1.1)
CurlyBraces(x0 = 1 * curlyBracesOffset,
x1 = lastDayOfWinter * curlyBracesOffset,
y0 = -0.1 * max(tempData$RAIN), y1 = -0.1 * max(tempData$RAIN),
position = 1,
bracesSize = 0.15 * max(tempData$RAIN),
label = paste0(
precipitationSums$winter[precipitationSums$Site == site &
precipitationSums$year == selectedYear],
" mm (winter)"),
labelSize = 1, labelDist = 1, labelAdj = 0)
CurlyBraces(x0 = (lastDayOfSummer + 1) * curlyBracesOffset,
x1 = 365 * curlyBracesOffset,
# x1 coordinate must be higher than 365 to cover the rest of the axis
y0 = -0.1 * max(tempData$RAIN), y1 = -0.1 * max(tempData$RAIN),
position = 1,
bracesSize = 0.15 * max(tempData$RAIN))
par(mar = c(0.1, leftMargin, 1.1, 0.2))
barplot(tempData$ETr,
ylab = yLabs[2])
plot(1:nrow(tempData), tempData$WATp * 100,
xlab = "", xaxt = 'n',
ylab = yLabs[3],
type = "l", lwd = 2, ylim = c(0, soil.FC * 110))
abline(h = soil.FC * 100, lty = 2)
abline(h = soil.WP * 100, lty = 2)
if (site == levels(factor(weather$Site))[1])
{
text(x = 5, y = 17.5, labels = "field capacity",
cex = fontScale * 0.8, font = 3, adj = 0)
text(x = 5, y = 4.5, labels = "wilting point",
cex = fontScale * 0.8, font = 3, adj = 0)
}
par(mar = c(3.5, leftMargin, 1.1, 0.2))
plot(1:nrow(tempData), tempData$ARID,
xlab = "day\n", ylab = yLabs[4],
type = "l", lwd = 2, ylim = c(0, 1))
}
rm(tempData)
dev.off()
## png
## 2
knitr::include_graphics(plotName)
Load output generated with NetLogo implementation of the SIMPLE crop model, which integrates the Weather and Soil Water Balance models:
simExample <- read.csv("Fig6_input/SIMPLE-crop-model experiment-table.csv", skip = 6)
wheatHarvestDay = simExample$"item.0.sugHarvestingDay"[1]
riceHarvestDay = simExample$"item.1.sugHarvestingDay"[1]
simExample_pars <- simExample[1, c(
"solar_annual.min", "solar_annual.max", "solar_daily.mean.fluctuation",
"temperature_annual.min.at.2m", "temperature_annual.max.at.2m",
"temperature_daily.mean.fluctuation",
"temperature_daily.lower.deviation", "temperature_daily.upper.deviation",
"precipitation_yearly.mean", "precipitation_yearly.sd",
"precipitation_daily.cum_plateau.value_yearly.mean",
"precipitation_daily.cum_inflection1_yearly.mean", "precipitation_daily.cum_inflection1_yearly.sd",
"precipitation_daily.cum_rate1_yearly.mean", "precipitation_daily.cum_rate1_yearly.sd",
"precipitation_daily.cum_plateau.value_yearly.sd", "precipitation_daily.cum_max.sample.size",
"precipitation_daily.cum_inflection2_yearly.mean", "precipitation_daily.cum_inflection2_yearly.sd",
"precipitation_daily.cum_rate2_yearly.mean", "precipitation_daily.cum_rate2_yearly.sd", "precipitation_daily.cum_n.sample"
)]
simExample <- simExample[, c(
"X.step.", "T", "RAIN", "solarRadiation", "ET_0",
"mean..WATp..of.patches", "mean..ARID..of.patches",
"mean..biomass..of.patches.with..position.crop.typesOfCrops...0.",
"mean..biomass..of.patches.with..position.crop.typesOfCrops...1.",
"mean..yield..of.patches.with..position.crop.typesOfCrops...0.",
"mean..yield..of.patches.with..position.crop.typesOfCrops...1."
)]
names(simExample) <- c(
"day", "T", "RAIN", "solarRadiation", "ETr",
"WATp", "ARID",
"wheat_biomass", "rice_biomass", "wheat_yield", "rice_yield"
)
# create vector with harvest yields per crop and year
years <- 1:floor(nrow(simExample) / 365) - 1
cropYields <- data.frame(
wheatHarvestDays = (years * 365 + wheatHarvestDay),
riceHarvestDays = (years * 365 + riceHarvestDay),
wheat = unique(simExample$wheat_yield),
rice = unique(simExample$rice_yield)[-1])
cropYields$wheat[1] <- NA # ignore first zero yield in wheat
Show table with parameter values:
parvalues <- c()
for (i in 1:length(simExample_pars))
{
parvalues <- c(parvalues, simExample_pars[[i]])
}
parvalues <- cbind(names(simExample_pars), parvalues)
knitr::kable(parvalues,
format = "html",
col.names = c("parameter", "values"),
align = c("l", "r"))
| parameter | values |
|---|---|
| solar_annual.min | 3 |
| solar_annual.max | 7 |
| solar_daily.mean.fluctuation | 1 |
| temperature_annual.min.at.2m | 15 |
| temperature_annual.max.at.2m | 40 |
| temperature_daily.mean.fluctuation | 5 |
| temperature_daily.lower.deviation | 5 |
| temperature_daily.upper.deviation | 5 |
| precipitation_yearly.mean | 400 |
| precipitation_yearly.sd | 130 |
| precipitation_daily.cum_plateau.value_yearly.mean | 0.1 |
| precipitation_daily.cum_inflection1_yearly.mean | 40 |
| precipitation_daily.cum_inflection1_yearly.sd | 20 |
| precipitation_daily.cum_rate1_yearly.mean | 0.15 |
| precipitation_daily.cum_rate1_yearly.sd | 0.02 |
| precipitation_daily.cum_plateau.value_yearly.sd | 0.05 |
| precipitation_daily.cum_max.sample.size | 10 |
| precipitation_daily.cum_inflection2_yearly.mean | 200 |
| precipitation_daily.cum_inflection2_yearly.sd | 20 |
| precipitation_daily.cum_rate2_yearly.mean | 0.05 |
| precipitation_daily.cum_rate2_yearly.sd | 0.01 |
| precipitation_daily.cum_n.sample | 200 |
Load function for marking the end of each year:
markEndYears <- function(lengthOfData, offset = 1){
for (i in 1:lengthOfData)
{
if (i %% (365 * offset) == 0)
{
abline(v = i, lty = 3)
}
}
}
Plot all selected variables:
#--------
# To create a png file:
plotName = "Fig6.png"
fontScale = c(2, 2, 1.2, # c(cex.lab, legend, cex.axis)
3, 2) # side annotation name and "model"
png(plotName, width = 1000, height = 1200)
#---------
# alternatively, to create eps file:
# plotName = "Fig6.eps"
#
# fontScale = c(1.8, 1.8, 1.2, # c(cex.lab, legend, cex.axis)
# 2.8, 1.8) # side annotation name and "model"
#
# extrafont::loadfonts(device = "postscript")
# grDevices::postscript(file = plotName,
# pointsize = 10,
# width = 10,
# height = 12,
# horizontal = FALSE,
# paper = "special",
# onefile = FALSE,
# family = "sans",
# colormodel = "cmyk")
#---------
layout(matrix(c(1:8, rep(9, 4), rep(10, 2), 11, 12),
nrow = 8, ncol = 2),
heights = c(1, 1, 1, 1, 1, 1, 1.3, 0.26),
widths = c(1, 0.15))
yLabs <- c(expression(paste(" Solar\nRadiation (", MJ/m^-2, ")")),
"Temperature (ºC)",
"Precipitation (mm)",
"ETr (mm)", "WATp (%)", "ARID index")
leftMargin <- 7
# 1. Solar radiation
par(mar = c(0.1, leftMargin, 1.1, 0.2),
cex.lab = fontScale[1],
cex.axis = fontScale[3])
plot(1:nrow(simExample), simExample$solarRadiation, type = "l",
xlab = "", xaxt = 'n',
ylab = yLabs[1])
markEndYears(nrow(simExample))
# 2. Temperature
plot(1:nrow(simExample), simExample$"T", type = "l",
xlab = "", xaxt = 'n',
ylab = yLabs[2])
#lines(1:nrow(simExample), simExample$)
markEndYears(nrow(simExample))
# 3. Precipitation
par(mar = c(1.1, leftMargin, 2.1, 0.2))
barplot(simExample$RAIN,
ylab = yLabs[3])
markEndYears(nrow(simExample), offset = 1.2) ; abline(v = 2190 * 1.2, lty = 3)
# not sure why, but barplot() x coordinates do not behave as in plot()
# 4. Reference evapotranspiration
par(mar = c(0.1, leftMargin, 1.1, 0.2))
barplot(simExample$ETr,
ylab = yLabs[4])
markEndYears(nrow(simExample), offset = 1.2) ; abline(v = 2190 * 1.2, lty = 3)
# not sure why, but barplot() x coordinates do not behave as in plot()
# 5. Soil water content (%)
plot(1:nrow(simExample), simExample$WATp * 100,
xlab = "", xaxt = 'n',
ylab = yLabs[5],
type = "l", lwd = 2, ylim = c(0, soil.FC * 110))
abline(h = soil.FC * 100, lty = 2)
abline(h = soil.WP * 100, lty = 2)
markEndYears(nrow(simExample))
text(x = 2000, y = 17.5, labels = "field capacity",
cex = fontScale[1], font = 3)
text(x = 2000, y = 4.5, labels = "wilting point",
cex = fontScale[1], font = 3)
# 6. ARID index
plot(1:nrow(simExample), simExample$ARID,
xlab = "", xaxt = 'n',
ylab = yLabs[6],
type = "l", lwd = 2, ylim = c(0, 1))
markEndYears(nrow(simExample))
# 7. Crop biomass and yield
par(mar = c(3.5, leftMargin, 1.1, 0.2))
plot(1:nrow(simExample), simExample$wheat_biomass, type = "l",
xlab = "", col = "darkred", lwd = 2,
ylab = expression(paste("Biomass (", g/m^-2, ")")))
lines(1:nrow(simExample), simExample$rice_biomass, col = "darkblue", lwd = 2)
lines(cropYields$wheatHarvestDays, cropYields$wheat, col = "darkred", lwd = 2, lty = 3, pch = 4, cex = 3, type = "o")
lines(cropYields$riceHarvestDays, cropYields$rice, col = "darkblue", lwd = 2, lty = 3, pch = 4, cex = 3, type = "o")
markEndYears(nrow(simExample))
text(x = 2050, y = 700, labels = "harvested",
cex = fontScale[1], font = 3)
arrows(x0 = 2050, x1 = 1950, y0 = 620, y1 = 450, length = 0.1, lwd = 2)
mtext("day", side = 1, line = 3, cex = fontScale[3])
CurlyBraces(x0 = 0.02, x1 = 364.08,
y0 = -100, y1 = -100,
position = 1,
bracesSize = 100)
mtext("year (365 days)", side = 1, line = 3,
cex = fontScale[3], adj = 0.06, font = 3)
# 8. Legend (crops)
par(mar = c(0,0,0,0), cex = fontScale[2])
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
legend(x = 0.05, y = 1.1,
legend = c('Triticum aestivum (Yecora Rojo)', 'Oryza sativa (IR72)'),
col = c('darkred', 'darkblue'), lty = c(1, 1), lwd = 2,
horiz = T, bty = "n")
# 9. Weather model bracket
par(mar = c(0,0,0,2), cex = fontScale[2])
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
CurlyBraces(x0 = 0.01, x1 = 0.01,
y0 = -0.02, y1 = 1.01,
position = 4,
bracesSize = 0.66,
label = "Weather\n",
labelSize = fontScale[4])
mtext(text = "model",
side = 4, line = 0.33, cex = fontScale[5])
# 10. Soil Water Balance model bracket
par(mar = c(0,0,0,2), cex = fontScale[2])
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
CurlyBraces(x0 = 0.01, x1 = 0.01,
y0 = -0.02, y1 = 1,
position = 4,
bracesSize = 0.66,
label = "Soil Water Balance\n",
labelSize = fontScale[4])
mtext(text = "model",
side = 4, line = 0.33, cex = fontScale[5])
# 11. Crop model bracket
par(mar = c(0,0,0,2), cex = fontScale[2])
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
CurlyBraces(x0 = 0.01, x1 = 0.01,
y0 = 0.1, y1 = 1,
position = 4,
bracesSize = 0.66,
label = "Crop\n",
labelSize = fontScale[4])
mtext(text = "model",
side = 4, line = 0.33, cex = fontScale[5])
# 12. empty
plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
dev.off()
## png
## 2
knitr::include_graphics(plotName)
Load packages:
library(grid)
library(gridExtra)
## Warning: package 'gridExtra' was built under R version 3.6.3
library(png)
Defining terrain random number generator seed (must be available inside “terrains” folder:
listOfSeeds <- c(32, 298, 304, 499)
Setting dimensions:
figScale = 2
terrainDim = 450 * figScale
topHeight = 100 * figScale
bottomHeight = 150 * figScale
leftWidth = 100 * figScale
internalMargin = 5 * figScale
figWidth = terrainDim * 3 + leftWidth + 2 * internalMargin
figHeight = terrainDim * 4 + topHeight + bottomHeight + 3 * internalMargin
Create margin header images:
headers <- c("seed",
"Elevation",
"Soil texture type",
"Ecological community\n(cover type)")
namesInFile <- c("seed", "elev", "soil", "ecol")
for (i in 1:length(headers))
{
if (i == 1) { width = leftWidth } else { width = terrainDim }
png(paste0("Fig7_input/Fig7_", namesInFile[i], ".PNG"),
width = width,
height = topHeight)
par(mar = rep(0, 4))
plot(0:1, 0:1, type = "n", axes = FALSE, ann = FALSE)
text(0.5, 0.5, cex = figScale * 3.7,
labels = headers[i])
dev.off()
}
for (seed in listOfSeeds)
{
png(paste0("Fig7_input/Fig7_seed=", seed, ".PNG"),
width = leftWidth,
height = terrainDim)
par(mar = rep(0, 4))
plot(0:1, 0:1, type = "n", axes = FALSE, ann = FALSE)
text(0.5, 0.5, cex = figScale * 5, adj = .5,
labels = seed)
dev.off()
}
Create legend images:
# code converted from NetLogo (maximum = 255)
elev <- 100/255 + (155/255) * seq(0, 1, length.out = 6)
png("Fig7_input/Fig7_legend_elev.PNG",
width = terrainDim,
height = bottomHeight)
par(mar = c(2,2,1,2), cex = figScale * 2)
barplot(rep(100, 6),
col = rgb(elev - (100/255), elev, 0),
axes = FALSE)
mtext("min.", side = 1, line = 0.1, adj = 0, cex = figScale * 2)
mtext("max.", side = 1, line = 0.1, adj = 1, cex = figScale * 2)
#mtext("elevation", side = 1, line = 1.2, adj = .5, cex = figScale * 3)
dev.off()
soilTextureTypes <- c("Sandy loam", "Loam",
"Sandy clay loam", "Clay loam",
"Sandy clay", "Clay")
soilTextureTypes_colors <- c("#9d6e48", "#eded31",
"#a71b6a", "#345da9",
"#7c50a4", "#2d8dbe")
xy <- list(c(-0.1, 1), c(0.55, 1),
c(-0.1, 0.7), c(0.55, 0.7),
c(-0.1, 0.4), c(0.55, 0.4))
png("Fig7_input/Fig7_legend_soil.PNG",
width = terrainDim,
height = bottomHeight)
par(mar = rep(0, 4))
plot(0:1, 0:1, type = "n", axes = FALSE, ann = FALSE)
for (i in 1:length(soilTextureTypes))
{
legend(xy[[i]][1], xy[[i]][2],
legend = soilTextureTypes[i],
fill = soilTextureTypes_colors[i],
bty = "n", adj = 0.1,
cex = figScale * 2.5, pt.cex = figScale * 5)
}
dev.off()
ecologicalCommunityTypes <- c("grassland", "wood-grass",
"shrubland", "woodland")
ecologicalCommunityTypes_colors <- c("#f16a15", "#eded31",
"#59b03c", "#1d9f78")
xy <- list(c(0, 0.95), c(0.45, 0.95),
c(0, 0.5), c(0.45, 0.5))
png("Fig7_input/Fig7_legend_ecol.PNG",
width = terrainDim,
height = bottomHeight)
par(mar = rep(0, 4))
plot(0:1, 0:1, type = "n", axes = FALSE, ann = FALSE)
for (i in 1:length(ecologicalCommunityTypes))
{
legend(xy[[i]][1], xy[[i]][2],
legend = ecologicalCommunityTypes[i],
fill = ecologicalCommunityTypes_colors[i],
bty = "n", adj = 0.1,
cex = figScale * 2.5, pt.cex = figScale * 5)
}
dev.off()
Create blank image:
png("Fig7_input/Fig7_blank.PNG",
width = 10,
height = 10)
## Warning in png("Fig7_input/Fig7_blank.PNG", width = 10, height = 10): 'width=10,
## height=10' are unlikely values in pixels
dev.off()
## png
## 2
Load input images:
elevTerrains <- paste0("Fig7_input/terrains/terrain_random_w=50_h=50_a=C#_fill-sinks=true_seed=",
listOfSeeds ,".PNG")
soilTerrains <- paste0("Fig7_input/terrains/terrain_random_w=50_h=50_a=C#_fill-sinks=true_seed=",
listOfSeeds ,"_soilTextureTypes.PNG")
ecolTerrains <- paste0("Fig7_input/terrains/terrain_random_w=50_h=50_a=C#_fill-sinks=true_seed=",
listOfSeeds ,"_coverTypes.PNG")
# it is important that these plots were saved as "PNG not "png"
listOfPlots <- c("Fig7_input/Fig7_seed.PNG",
"Fig7_input/Fig7_elev.PNG",
"Fig7_input/Fig7_blank.PNG",
"Fig7_input/Fig7_soil.PNG",
"Fig7_input/Fig7_blank.PNG",
"Fig7_input/Fig7_ecol.PNG",
# ---
paste0("Fig7_input/Fig7_seed=", listOfSeeds[1], ".PNG"),
elevTerrains[1],
"Fig7_input/Fig7_blank.PNG",
soilTerrains[1],
"Fig7_input/Fig7_blank.PNG",
ecolTerrains[1],
# ---
rep("Fig7_input/Fig7_blank.PNG", 6),
# ---
paste0("Fig7_input/Fig7_seed=", listOfSeeds[2], ".PNG"),
elevTerrains[2],
"Fig7_input/Fig7_blank.PNG",
soilTerrains[2],
"Fig7_input/Fig7_blank.PNG",
ecolTerrains[2],
# ---
rep("Fig7_input/Fig7_blank.PNG", 6),
# ---
paste0("Fig7_input/Fig7_seed=", listOfSeeds[3], ".PNG"),
elevTerrains[3],
"Fig7_input/Fig7_blank.PNG",
soilTerrains[3],
"Fig7_input/Fig7_blank.PNG",
ecolTerrains[3],
# ---
rep("Fig7_input/Fig7_blank.PNG", 6),
# ---
paste0("Fig7_input/Fig7_seed=", listOfSeeds[4], ".PNG"),
elevTerrains[4],
"Fig7_input/Fig7_blank.PNG",
soilTerrains[4],
"Fig7_input/Fig7_blank.PNG",
ecolTerrains[4],
# ---
"Fig7_input/Fig7_blank.PNG",
"Fig7_input/Fig7_legend_elev.PNG",
"Fig7_input/Fig7_blank.PNG",
"Fig7_input/Fig7_legend_soil.PNG",
"Fig7_input/Fig7_blank.PNG",
"Fig7_input/Fig7_legend_ecol.PNG"
)
grobs <- lapply(lapply(listOfPlots, png::readPNG), grid::rasterGrob)
Create Figure 7:
widths = 50 * c(leftWidth/figWidth,
terrainDim/figWidth,
internalMargin/figWidth,
terrainDim/figWidth,
internalMargin/figWidth,
terrainDim/figWidth) * figScale
heights = (figHeight/figWidth) * 50 * c(topHeight/figHeight,
terrainDim/figHeight,
internalMargin/figHeight,
terrainDim/figHeight,
internalMargin/figHeight,
terrainDim/figHeight,
internalMargin/figHeight,
terrainDim/figHeight,
bottomHeight/figHeight) * figScale
lay <- rbind(1:6, 7:12, 13:18, 19:24, 25:30, 31:36, 37:42, 43:48, 49:54)
png("Fig7.png", width = figWidth, height = figHeight)
gridExtra::grid.arrange(grobs = grobs,
layout_matrix = lay,
widths = unit(widths, "cm"),
heights = unit(heights, "cm"))
dev.off()
## png
## 2
knitr::include_graphics("Fig7.png")
terrainsData <- read.csv(
paste0("Fig7_input/terrains/terrain_random_w=50_h=50_a=C#_fill-sinks=true_seed=", listOfSeeds[1] ,".csv"), skip = 8, nrows = 1)
for (i in 2:length(listOfSeeds))
{
terrainsData <- rbind(terrainsData,
read.csv(
paste0("Fig7_input/terrains/terrain_random_w=50_h=50_a=C#_fill-sinks=true_seed=", listOfSeeds[i] ,".csv"), skip = 8, nrows = 1))
}
terrainsData <- as.data.frame(t(terrainsData[, c(89, 18:39,
94, 96, 99:102, 104:107, 111,
11:16)]))
#names(terrainsData) <- listOfSeeds
Show table with parameter values:
knitr::kable(terrainsData,
format = "html",
col.names = rep("", 4),
align = c("l", "c", "c", "c"))
| randomseed | 32 | 298 | 304 | 499 |
| elev_algorithm.style | “C#” | “C#” | “C#” | “C#” |
| elev_featureanglerange | 10.355587 | 13.886082 | 1.513911 | 15.282367 |
| elev_inversioniterations | 1.32550504 | 0.61290965 | 0.96763471 | 0.05134417 |
| elev_noise | 0.4589207 | 4.4218771 | 0.9439117 | 0.7888913 |
| elev_numdepressions | 2 | 9 | 2 | 5 |
| elev_numprotuberances | 2 | 1 | 8 | 8 |
| elev_numranges | 16 | 52 | 69 | 7 |
| elev_numrifts | 67 | 48 | 61 | 71 |
| elev_rangeaggregation | 0.8162051 | 0.2617838 | 0.4092565 | 0.1425585 |
| elev_rangeheight | 36.833480 | 5.124392 | 17.353900 | 30.715513 |
| elev_rangelength | 3311 | 769 | 1418 | 2718 |
| elev_riftaggregation | 0.59655344 | 0.04929822 | 0.18191958 | 0.73627220 |
| elev_riftheight | -30.45546 | -22.91505 | -43.80874 | -32.58051 |
| elev_riftlength | 3217 | 3021 | 3085 | 2546 |
| elev_smoothingradius | 3.464823 | 3.464823 | 3.464823 | 3.464823 |
| elev_smoothstep | 1 | 1 | 1 | 1 |
| elev_valleyaxisinclination | 0.7037125 | 0.7573202 | 0.0565110 | 0.6901128 |
| elev_valleyslope | 0.011788025 | 0.002599597 | 0.004806264 | 0.010067643 |
| elev_xslope | 0.004417135 | 0.005646562 | 0.007325109 | 0.007044657 |
| elev_yslope | 0.005514878 | 0.006371659 | 0.006696862 | 0.009834952 |
| flow_do.fill.sinks | true | true | true | true |
| flow_riveraccumulationatstart | 34790950 | 18916696 | 20365916 | 43365595 |
| soil_depthnoise | 4.645497 | 40.531347 | 38.775613 | 38.513392 |
| soil_formativeerosionrate | 1.685376 | 2.807310 | 1.731212 | 1.152460 |
| soil_max.clay | 75.65757 | 61.35114 | 64.55615 | 90.23157 |
| soil_max.sand | 94.74552 | 99.87023 | 97.47498 | 96.08472 |
| soil_max.silt | 99.96852 | 49.18063 | 64.34818 | 95.20315 |
| soil_maxdepth | 514.8621 | 230.2860 | 387.0962 | 297.2467 |
| soil_min.clay | 16.24237 | 61.34954 | 38.91594 | 73.28332 |
| soil_min.sand | 90.25214 | 99.54846 | 54.11990 | 88.72363 |
| soil_min.silt | 45.20185 | 31.31848 | 59.64769 | 91.84854 |
| soil_mindepth | 241.5896 | 194.5769 | 149.2301 | 147.3834 |
| soil_texturenoise | 1.606241 | 6.412468 | 1.963127 | 4.493502 |
| ecol_brushfrequencyinflection | 20.10049 | 14.47833 | 12.84533 | 26.76113 |
| ecol_brushfrequencyrate | 0.01765001 | 0.09451887 | 0.03295388 | 0.09052492 |
| ecol_grassfrequencyinflection | 4.0538838 | 3.2674253 | 1.9392423 | 0.5540513 |
| ecol_grassfrequencyrate | 0.008029435 | 0.133515886 | 0.047298747 | 0.109743989 |
| ecol_woodfrequencyinflection | 41.68921 | 54.95257 | 30.47011 | 32.99336 |
| ecol_woodfrequencyrate | 0.005506242 | 0.008327215 | 0.071837624 | 0.066332965 |